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Abstract 

The design of general finite multi-server queueing networks is a challenging problem that arises in many real-life situations, 
including computer networks, manufacturing systems, and telecommunication networks. In this paper, we examine the 
optimal routing problem in arbitrary configured acyclic queueing networks. The performance of the finite queueing 
network is evaluated with a known approximate performance evaluation method and the optimization is done by means of 
a heuristics based on the Powell algorithm. The proposed methodology is then applied to determine the optimal routing 
probability vector that maximizes the throughput of the queueing network. We show numerical results for some networks 
to quantify the quality of the routing vector approximations obtained. 



Citation: van Woensel T, Cruz FRB (2014) Optimal Routing in General Finite Multi-Server Queueing Networks. PLoS ONE 9(7): e102075. doi:10.1371/joumal.pone. 
0102075 

Editor: Gerardo Adesso, University of Nottingham, United Kingdom 
Received April 15, 2014; Accepted June 13, 2014; Published July 10, 2014 

Copyright: © 2014 van Woensel, Cruz. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits 
unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 

Data Availability: The authors confirm that all data underlying the findings are fully available without restriction. All relevant data are within the paper and its 
Supporting Information files. 

Funding: FRBC was supported by the Conselho Nacional de Desenvolvimento Cientffico e Tecnologico (CNPq) of the Ministry for Science and Technology of 
Brazil, grant number 303388/2010-2, and Fundacao de Amparo a Pesquisa do Estado de Minas Gerais (FAPEMIG), grant number CEX-PPM-00071-12. The funders 
had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. 

Competing Interests: The authors have declared that no competing interests exist. 

* Email: T.v.Woensel@tue.nl 



Introduction 

The design of networks with random arrivals, random service 
times, multiple servers, and a finite number of buffer spaces is a 
challenging problem that arises in many real-life situations, e.g. in 
manufacturing systems [1,2], computer networks [3,4], public 
services [5,6], call centers [7,8], pedestrian and vehicular traffic 
[9-12], among other situations. This problem is challenging 
because finite queueing networks are notoriously difficult to 
analyze analytically, and closed form expressions are not easily 
constructed for the performance measures of such systems. Also 
note that the underlying network design problems involved are 
usually very hard to solve. 

In fact, there are several distinct network design optimization 
problems associated with finite queueing networks. According to 
Daskalaki & Smith [13] the optimal network design problem can 
be split up into three interrelated optimization problems: 

1 . The optimal topology problem (OTOP), which deals with 
decisions of the design of the network itself, that is, the number 
of nodes (e.g. workstations, warehouses, delivery points, etc.) 
and arcs (e.g. corridors, conveyors, escalators, etc.) and the 
general configuration of these two components; 

2. The optimal routing problem (OROP), which deals with 
determining the routing probabilities in the network defined by 
the first problem; 

3. The optimal resource allocation problem (ORAP), which deals 
with the optimal allocation of the scarce resources in the 
network, e.g. the number of buffers (i.e., the buffer allocation 



problem, BAP) and the number of servers (i.e., the server 
allocation problem, CAP). 

These three problems are challenging and difficult optimization 
problems. For an arbitrary topology, the OTOP is shown to be 
MV — hard [14], and the same is conjectured for the general 
ORAP [15]. 

Previous work focused mainly on the ORAP in open finite 
acyclic queueing network settings. Both BAP and CAP are 
probably among the most well-known optimal resource allocation 
problems [16]. For instance, Cruz et al. [17] and Smith et al. [18] 
looked into the BAP, both in a single and in a multi-server setting, 
and Smith et al. [19] proposed algorithms to solve the CAP. 
However, the routing probabilities are usually assumed to be 
known beforehand for BAP and CAP [20]. 

The overall research objective of this paper is to build a relevant 
model and solution methodology for the system's throughput 
maximization problem. In this paper, we focus on optimizing the 
routing probabilities through the queueing network, i.e. the 
OROP. A similar research question is tackled by Daskalaki & 
Smith [13] in which they evaluated the joint effect of buffer 
allocation and routing on the throughput. Earlier, Gosavi & Smith 
[21] focused on the routing optimization problem related to the 
overall objective of throughput maximization. The common 
ground of both papers is that they used queueing networks with 
single servers and exponential service times [13,21]. Kerbache & 
Smith [22] considered, for different product classes, the optimal 
routes conjoint with a variant of the optimal topology problem to 
determine the connected arcs in a single server queueing network 
setting. Distinct from Daskalaki & Smith [13] and Gosavi & Smith 
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algorithm 

/* Step 1: Initialization */ 

1.1 re&dV(V,A) 

/* initialize the routing probability vector */ 

1.2 fc^O 

1.3 «W=«g*\v(t,.7)6i4 

/* evaluate the initial solution with the GEM */ 

1.4 8« «- 9(aW) 

/* Step 2: Optimization & Performance Evaluation */ 
/* generate new solution using Powell */ 

2.1 k<-Jfe + l 

2.2 d$ <- PoweU(a^" 1) ,e( fc - 1 )) ! V(i,j") £ A 

/* evaluate the current solution with the GEM */ 

2.3 ew <- e(a<*>) 

2.4 if |9 W - e^- 1 ^ > £ then goto 2.1 
/* Step 3: Print Results */ 

3.1 print a< fc > and 9^ 
end algorithm 

Figure 1. Structured overview of the methodology. 

doi:10.1371/journal.pone.0102075.g001 

[21] is that Kerbache & Smith [22] considered general arrival 
times, general service times, and single server queues. Secondly, 
Gosavi & Smith [21] did not consider the general expansion 
method (GEM) in their analysis as the evaluation tool [23]. 

Specifically, we examine the OROP, by means of a combina- 
tion of the GEM and a heuristic based on the Powell algorithm 
[24], specifically for acyclic networks of M/ G/c/K queues, which 
in Kendall notation means a queueing system with Markovian 
arrival rates, General service times, c parallel servers, and a total 
capacity of K items, including those items in service (practical 
applications to M/G/c/K queueing networks include manufac- 
turing and service systems [25] and transportation and material 
handling systems [26]). The results are compared to simulations to 
attest for the quality of the routing vectors obtained. Besides, 
another important contribution of this paper is to present helpful 
approximations to swift managerial decisions regarding the 
optimal routing probability vectors to maximize the overall 
throughput in a network of finite general-service queues. We also 
present important empirical arguments to show that these 
approximations are effective. 

This paper is organized as follows. The next section describes in 
detail the mathematical model formulation considered and 
elaborates on both the performance evaluation tool and on the 
optimization procedure. In the following section detailed results 
are given for the problem on hand. Finally, the last section 
concludes this paper with final remarks and topics for future work 
in the area. 

Materials and Methods 

Mathematical Programing Formulation 

The problem is defined on a digraph T> = (V ,A), in which V is 
the set of vertexes (finite queues) and A is the set of arc 
(connections between the queues). Each vertex (queue) is 
characterized by Poisson arrivals, general service, and multiple 
servers. Mathematically, the optimal routing problem can be 
formulated as follows. 

(OROP): 

max 0(a), (1) 



subject to: 

0<a, v <l, V(i,j)eA, (2) 
£ay = l, Vi ' eK > ( 3 ) 

Y/e*") 

in which 0(a) is the throughput, which is the objective that must 
be maximized, a the optimal routing probability matrix, a = [a j\ , 
i.e. the matrix that maximizes the objective function of this 
predefined network, and 5(f) is the set of succeeding vertexes of 
vertex i, that is, S(i) = {j\ (i, j) e A}. 

The throughput will thus be affected by the effective routing of 
jobs through the network, the variability of the service distribution, 
the number of servers, and the number of buffers. It should be 
clear that the above model is a highly difficult stochastic 
programming problem to handle due to the inherent non-linearity 
of the objective function combined with the absence of any closed- 
form expressions for the throughput 0(a). 

Proposed Algorithm 

Figure 1 presents the proposed algorithm to solve the OROP in 
order to provide more insights into the interaction between the 
performance evaluation tool and the optimization method. 

The initial routing probability vector is conveniently set to the 
inverse of the number of nodes after a split, 

«<f> =ivy^, (4) 

in which is the number of succeeding nodes to node i, that is, 
the cardinality of set A(i). The optimization step itself is an 
iteration in which new solutions are generated following the Powell 
logic until convergence, that is, until the difference in 0, 
A0 = (0 </l) — is less than a predefined tolerance e. 

The Powell algorithm can be described as an unconstrained 
optimization procedure that does not require the calculation of 
first derivatives of the function. Numerical examples have shown 
that the method is capable of minimizing a function with up to 
twenty variables [24,27]. The Powell method locates the minimum 
of a non-linear function /(x) by successive uni-dimensional 
searches from an initial starting point x' ' along a set of conjugate 
directions. These conjugate directions are generated within the 
procedure itself. The Powell method is based on the idea that if a 
minimum of a non-linear function f(x) is found along p conjugate 
directions in a stage of the search, and an appropriate step is made 
in each direction, the overall step from the beginning to the p-th 
step is conjugate to all of the p sub-directions of the search. We 
have seen remarkable success in the past with coupling Powell 
algorithm and the GEM [19], We discuss the GEM in detail now, 
which is also the method we used to obtain the performance 
measures for the problem studied in this paper. 

Performance Evaluation 

In previous papers (see e.g. Kerbache & Smith [23,28]) the GEM 
has been successfully used to evaluate the performance measures of 
acyclic networks of finite queues. The GEM is a robust and effective 
approximation technique that is basically a combination of repeated 
trials and node-by-node decomposition in which each queue is 
analyzed separately and then corrections are made in order to take 
into account the interrelation between the queues in the network. 
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Figure 2. The generalized expansion method. 

doi:1 0.1 371 /journal.pone.01 02075.g002 



The GEM has three stages, Network Reconfiguration, Parameter 
Estimation, and Feedback Elimination, to be described as follows. 

Stage I: Network Reconfiguration. The first step in the 
GEM involves reconfiguring the network. An artificial vertex hj is 
added preceding each finite vertex j in the network. The artificial 
vertex is added to register the blocked customers at node j and is 
modeled as an M / Gj oo queue, as shown in Figure 2. 

When an entity leaves vertex i, vertex / may be blocked with 
probability Pk p or unblocked, with probability (1— Pk,)- Under 
blocking, the entities are rerouted to vertex hj for a delay while 
node j is busy. Vertex hj helps to accumulate the time an entity has 
to wait before entering vertex j and to compute the effective arrival 
rate to vertex /. In other words, the GEM ultimate goal is to 
provide an approximation scheme to update the service rates at the 
upstream vertex i to take into account all blocking after service 
caused by the downstream vertex j : 



(5) 



in which p t is the exponential service rate at vertex i, p^. is the 
blocking probability of finite queue j of size Kj, p},. is the corrected 
exponential service rate at the artificial vertex hj, and Jl i is the 
updated service rate at vertex i. As a final note, an important point 
about this process is that we do not physically modify the networks, 
only represent the expansion process through the software. 

Stage II: Parameter Estimation. In the second stage, the 
parameters px, pk, an d Pi, must be estimated, which is done 
essentially by utilizing known results for queueing theory. Index j is 
omitted for simplicity. 

Pk : Analytical results from the M /M/c/K queue provide 
the following expression for the blocking probability p k ■ 



Px = 



1 

c K ~ c c\ 



Po, 



(6) 



in which 



PO = 



c-1 



2-<n\ \p) c\ l-A/(c/i) 



_«=o 
for \/(cp) # 1 



, (7) 



However, the interest is on M/G/c/K models, for 
which there is not exact closed form expression for p^ . 
Then approximations must be used and Kimura's [29] 
two moment approximation has proven to be very 
effective in similar cases [25,30]. For example, let us fix 
c = 2 and the following closed form expression can be 
developed from Equation (6), for the optimal buffer size 
Bm=K — 2 for Markovian M/M/2/K queues, as a 
function of the blocking probability: 



In 



Bm '- 



1 p K (2p + X) 
22p — X+p K X 



(8) 



The following Kimura's two moment approximation can 
be used to approximate the optimal buffer size B e {s 2 ) of 
a general service M/G/2/K queue: 



B e (s 2 )-- 



^m+NINT 



1 



yfpBi 



(9) 



in which s 2 is the squared coefficient of variation of the 
service time distribution at the queue, p = \/(cp) is the 
traffic intensity, Bm is the exact buffer size for a 
respective Markovian system, and NINTj x) is the 
nearest integer to x. Now, if we invert Equation (9) to 



«1,2 
M/G/4/20 



A, 




Figure 3. Basic split network B1. 

doi:1 0.1 371/joumal.pone.01 02075.g003 
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(a) effect of Ai versus qi,2 on throughput 0 
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(b) contour plot for Ai versus ai t 2 
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Figure 4. The shape of the objective function. 

doi:1 0.1 371 /journal, pone.01 02075.g004 

solve for pn we achieve: 



PK 



Pk - 



2p 



-2p 













-St} 


(2/i 


-A) 














A + 


2/i + A 



(10) 



This is a process that can be extended for c>2. In fact, 
expressions for px of up to c = 500 are available [30]. 
Another expressions, for c = 3, ...,10, are included in 
the software so that we have a complete set of blocking 
probabilities for c 6 {1, ... ,10}. 

Since there is no closed form solution for this quantity the 
following approximation is used given by Labetoulle & 
Pujolle [31] obtained using diffusion techniques. 



M/G/c 2 /2 



Oil, 2 



A i 




M/G/ci/20 / M/G/cq/2 
*f 



Figure 5. Basic split network B2. 

doi:1 0.1 371/joumal.pone.01 02075.g005 
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+ A((,f-rf)-(rf-'-rf-')) r 



/(/ ,((,f + l-,f+l)-(,f- f f)) / 
in which r\ and T2 are the roots to the polynomial 

\ — (\+li h + Hj)x+ n h x 2 = 0, 

with X = Xj — A/,(l — p'g), in which Xj and A/, are the 
actual arrival rates to the finite and artificial holding 
nodes respectively. Furthermore, the arrival rate to the 
finite node j is given by: 

Xj = Xi(l—pK) = Xi — Xh- 



N Ol 



fi/, : The delay distribution of a blocked customer at the 
holding node has the same distribution as the remaining 
service time of the customer being serviced at the node 
doing the blocking. Using renewal theory, one can show 
that the remaining service time distribution has the 
following rate fi/, : 



VO «- 
vo ro CTi 



2/i 



L 



-<tV ! 



(12) 



in which a is the service time variance given by 
Kleinrock [32]. Notice that if the service time distribu- 
tion at the finite queue doing the blocking is exponential 
with rate jij, then: 



that is, the service time at the artificial node is also 
exponentially distributed with rate flj. When the service 
time of the blocking node is not exponential, then ji h will 
be affected by a 2 . 
Stage III: Feedback Elimination. This stage is simply to 

eliminate the feedback loop, by recomputing the service time at 

vertex h^. The updated service rate is given by: 



*0 LT) 0\ (N 

a\ a* oo oo 



ri,=(\-p' K )fl h . 



Summary. Similar equations can be established with respect 
to each of the finite vertexes (queues). Ultimately, we have 
simultaneous non-linear equations in variables pn, pk, and [if,, 
along with auxiliary variables such as pLj and A,-. Solving these 
equations simultaneously, we can compute all performance 
measures of the network. 



Numerical Results and Discussion 

The software implementation is currently in Fortran 77. The 
compiler used was the DIGITAL Visual Fortran, Version 6, with 
the IMSL Fortran 90 MP Library version 3.0 for Microsoft 
Windows, to solve the nonlinear equations from the GEM. In our 
implementation, we set s= 10~ 1,000 , which proved to be effective 
based on the experiments. We first discuss the shape of the 
objective function. Secondly, we will give more insights for a 
number of split structures. We end the numerical results with some 
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Figure 6. Network structure CI. 

doi:1 0.1 371 /journal.pone.01 02075.g006 

complex network structures. Please bear in mind that the range of 
possible experiments is exponential itself, so we have determined a 
very selected, but representative sample to present. 

Shape of the Objective Function 

It is interesting to analyze the shape of the objective function for 
the optimization problem described earlier. The case discussed 
here is defined as follows. We have a three node network with a 
split into two branches, as seen in Figure 3. The general 
parameters for the base case are C\ = 4, K\ = 20, and c, = 2 and 
Ki = 2, for i = 2,3. The criteria to select those parameters is such 
that the number of servers C\ and the total capacity K\ of node 1 is 
larger than the others as to prevent it from becoming a bottleneck. 

We are particularly interested in the relationship between the 
overall throughput © = 0i+02> the routing probability aj 2, the 
arrival rate \\ , and the squared coefficient of variation of node 2, 
s\. Consequently, we set fl t = 2, for all nodes, and ^=$1 = 1. The 
sensitivity of these settings on the throughput is not analyzed now, 
but is amongst others the subject of study in the next sections. 
Next, we enumerate all possible combinations for Ai, «i 2, and s\, 
and then analytically obtain the corresponding throughput 0, 
which is shown in Figure 4 (always on the vertical axis), as a 
function of Ai, a 12, and s\- 

Figure 4-(a) clearly shows that the arrival rate is interacting with 
the routing probability. For low arrival rates, the network has low 
utilization. Consequently, different routing probabilities do not 
result in large changes in throughput 0. On the other hand, for 
large arrival rates, Ai>5, one clearly sees an optimal point in 
regard to the routing probability. Due to the symmetrical structure 
considered, a 50% split seems to be optimal here. Figure 4-{b) 
looks into the joint effect of changing the squared coefficient of 
variation, s\, together with the routing probability a 1,2- Again the 
inverted U-shape effect with a maximum around the 50% routing 
probability is visible. Next to this, it is clear that increasing the 
squared coefficient of variation from 0 to 2 reduces the overall 
throughput 0 but has a smaller impact on throughput than the 
routing probability. Based on this simple network structure, it is 
evident that the routing probabilities and the squared coefficient of 
variation affect the throughput to a large extent. Consequently, 



correctiy setting the routing matrix a leads to significant gains in 
terms of throughput. 

Basic Split Networks 

In this section, we analyze further the two-branch network from 
Figure 3 and include in our analysis the three-branch network seen 
in Figure 5. We are interested in assessing the influence of the 
number of servers c,, total capacities Kf, service rates fj, ( , and 
squared coefficient of variation of the service times sf, V7 £ V, in 
the model OROP, Equations (1) - (3). We choose to start with 
these two variants of a basic split structure as, from a routing 
allocation point of view, splits are the key building blocks in a 
generally configured network. The nodes after the splits are the 
ones of interest here. The first buffer K\ = 20 is larger than the 
others (Ki = 2, i = 2,3,4) as this will help to prevent the first queue 
of becoming a bottleneck node. The arrival rate A; is set equal to 
values {3,5,7}. 

Split with Two Branches. The top part of Table 1 gives the 
results for a two-branch split networks for unbalanced service rates 
ft and different squared coefficients of variation s\ and S3. In these 
cases the service rate of node 2 is made either relatively lower 
(/.( 2 = 1 versus ^3 =2), or equal (fi 2 =2 versus /.t 3 = 2), or higher 
(fi 2 = 3 versus ^3 =2) than the service rate of node 3. The base 
cases (sets Bid to B If) show that the routing probability is roughly 
equal to 0.5 when the nodes after the split are identical (that is, 
same number of servers, capacities, service rates, and squared 
coefficient of variation). Moreover, this results appears to be 
insensitive to changes in the squared coefficient of variation of 
both nodes after the split. Of course, the throughput 0 is affected 
(reduced) due to the changes (increase) in the variability. Secondly, 
changing the service rate of node 2, \i 2 (and keeping all other 
parameters equal to the base case settings), clearly shows that the 
fast nodes receive preference over the slow nodes. For example, in 
sets Bla to Blc (i.e., when node 2 is slower than node 3) a lower 
routing probability is set to node 2 (0.3334) than the one to node 3 
(0.6666). 

Rather than changing the squared coefficient of variation of 
both nodes after the split, we evaluated some unbalanced cases 
where only node 2 is affected by a different squared coefficient of 
variation, ^ = {0.0,0.5,1.0,1.5,2.0} (sets Blj to Blx, Table 1), 
combined with (^,^2,^3) = {(2,1,1),(2,2,2),(2,3,3)}. For these 
cases, we observe that the unbalance caused by the squared 
coefficient of variation only slightly changes the routing probability 
compared to sets with equal squared coefficients of variation (sets 
Bll, Blq, and Blv, Table 1). This is a confirmation of what we 
observed when evaluating the objective function earlier in the 
previous section. As we are now focusing on the small scale 
networks, this conclusion does not mean that the squared 
coefficient of variation has little effect in general. It is interesting 
to see that the throughput value is reducing as the squared 
coefficient of variation goes up although the routing probability is 
changing to protect more against the uncertainty in the second 
node. This is more prevalent in highly loaded systems. 

For the two-branch split networks, we evaluated a number of 
routing vectors around the optimal routing obtained. Table 2 
shows that the algorithm seems to have reached the optimal 
allocation for the routing probabilities into nodes 2 and 3 (set Ble, 
Table 1). Of course, one might argue that the optimization is 
rather easy due to the symmetric setting of the parameters. 
Therefore, we did the same analysis for the same parameter 
settings but with a network with unbalance in the service rates (set 
Bib, Table 1), also seen in Table 2. 

In conclusion, we observed that in previous results the 
optimization algorithm tries to balance out the flow taking into 
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account the differences (in service rates and squared coefficient of 
variation) among the two nodes after the split, which is intuitively 
logical as this strategy leads to the highest throughput. Moreover, 
the methodology seems to always find the optimal routing vector. 

Split with Three Branches. Let us now turn to the three- 
branch split network, Figure 5. It would be interesting to see to 
what extent the optimization algorithm balances the flow over the 
three nodes after the split and to what extent this is affected by the 
characteristics of the different nodes after the split. Table 3 shows a 
selected set of experiments done for this specific case. 

Table 3 shows that for the complete symmetric case, that is, set 
B2c, Table 3, again the routing probabilities are symmetric, i.e. 
a,j = 0.333,V (ij) e A. For the unbalanced cases in the squared 
coefficient of variation (sets B2a, B2b, B2d, and B2e, Table 3), it 
can be observed that the routing probability into the two identical 
nodes (a i > 3 and 1*1,4) are close to each other. For the remaining 
asymmetrical cases (sets B2f to B2o, Table 3), again the same 
conclusion holds. The faster (either in high number of servers or 
service rates) or more reliable (in terms of low squared coefficient 
of variation) are the nodes, more favored they are, resulting in high 
routing probabilities into these nodes. 



<^ <j* <ji 



»- r- o 



Complex Networks 

The simple networks discussed so far are interesting as they 
make it possible to show the behavior and logic of the optimization 
model in the presence of one split. In this section, we will evaluate 
some different complex topologies with regard to their routing 
probabilities. The first complex network considered is an extension 
of the two- and three-branch split networks, as depicted in 
Figure 6. 

Table 4 gives an overview of a selected set of experiments for 
the structure CI. The initial setting is again a balanced case, that 
is, c\=5, K\=20, H\ =2, c,=2, K,=2, ft = 2, sf = l, for 
! = 2,3,...,7 (set Cla, Table 4). Additional set of experiments 
involves unbalancing the service rates ft, the squared coefficients 
of variation sj, and the number of servers C;, for nodes 6 and 7. 
With these experiments, we evaluate whether and how the 
methodology takes the characteristics of the complete sub-network 
after the split into account in determining the optimal routing 
vector. 

We set up the experiments in such a way that either there are 
slow nodes (experiments Clc, Cle, Clg and Clh) or slow 
subsystems consisting of three connected nodes (experiments Clb, 
Cld, Clf, and Cli). Based on Table 4, we observe that in general 
the slower part of the network tends to receive less flow due to a 
lower routing probability into the relevant part. When after the 
first split in node 1 there is the choice to go to either the fast or 
slow subsystem, the faster subsystem is preferred. This is very clear 
in experiments Clb, Cld, Clf, and Cli, when the routing 
probability always favors the fastest downstream subsystem. 
However, if the last nodes are different (experiments Clc, Cle, 
Clg, and Clh), the conclusion is different. In all these 
experiments, the first split is just exactly half. The imbalance in 
the last nodes (i.e. nodes 4, 5, 6, and 7 are different), is completely 
absorbed in the routing probability at the immediately preceding 
nodes (i.e. nodes 2 and 3). Interestingly, this effect did not 
propagate upstream and did not affect the routing at the first split. 
Again, we see that the effect of the squared coefficient of variation 
on the routing probability is smaller compared to the number of 
servers or the service rates. 

The second network structure C2 has a more general structure 
than the other networks, as seen in Figure 7. Nodes 10 and 12 can 
act as a bottieneck node which might become overloaded 
depending upon the specific parameters. It is then interesting to 
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M/G/c 3 /K 3 M/G/a 7 /K 7 




M/G/cg/Kg 



Figure 7. Network structure C2. 

doi:10.1371/journal.pone.0102075.g007 



see how the routing probabilities are adapted to avoid or to reduce 
the workload in these bottleneck nodes. 

Table 5 gives the results for a selection of parameter settings for 
network structure C2. The table shows that when node 10 
becomes the bottleneck, routing the jobs into the direction of node 
10 is avoided by reducing the routing probability at node 3 (always 
around 0.4) and node 5 (ranging between 0.0409 to 0.4751). On 
the other hand, when the characteristics of node 1 0 are such that it 
is not a botdeneck, then the routing to nodes 5 and 6 is almost 50/ 
50. Secondly, it is clear that if the last node (node 12) becomes the 
bottleneck, only the throughput will be reduced. 

Approximations for the Routing Probabilities 

From a managerial point of view, it is interesting to have some 
good easy approximations that can be used to quickly set the 
routing probabilities. A number of possible approximations for the 
routing probabilities in the arc (ij) e A, after a split of vertex i 
into ii; vertexes, can be considered. 



JX). 



1 



„P). 



E tij 

VM0 



■VftT) e A, 



y(i,j) e A, 



E c 



e A, 



E CjHj 



,V(v) e A, 



(13) 



(14) 



(15) 



(16) 



in which S(i) is the set of succeeding vertexes of vertex i, that is, 
S(i) = {j\ (ij) e A}. Notice that «,- is the cardinality of set S(i). 

The first approximation, Equation (13), is simple but does not 
use any information from the W; vertexes after the split. This 
approximation only provides an equal spread of the throughput 
over the succeeding vertexes. It is expected that this approxima- 
tion works well when the nodes after the split are very similar in 
terms of service rate, number of servers, and so on. The other 



approximations, Equations (14), (15), and (16), do take more 
information into account. Equation (16) is believed to be the most 
general as it combines information in regard to the speed and the 
number of servers. On the other hand, no information about the 
squared coefficient of variation is included in none of the 
approximations . 

Tables 6 and 7 show that the performance of approximation acf 1 
improves as the network becomes more unbalanced (for instance, 
cases B 1 a-B 1 c are unbalanced, as defined in Table 1 , and for these 
cases the smallest A% found is for af\ A% = 0.00%, on the other 
hand cases B 1 d-B 1 f are balanced and for them all A% is equal to 
0.00%). This approximation of course takes into account the most 
information from the nodes after the split (not taking into account 
the squared coefficient of variation). If the nodes after the split are 
more alike (balanced) then the second approximation becomes 
favorable. On the other hand the first approximation a? is 
performing acceptable as well and could be preferred due to the 
easy implementation. 

Conclusions and Final Remarks 

In this paper, we examined the optimal routing problem in open 
finite acyclic queueing networks with a given general topology and 
multiple generally distributed servers. We determined the optimal 
routing probability vector that maximizes the throughput of an 
arbitrary configured network via a combination of the Generalized 
Expansion Method and Powell optimization tool. We presented 
numerical results showing the merits of the approach. Approxi- 
mations for the routing probability vector are also presented and 
evaluated. 

We have considered here only the throughput as the main 
performance measure. It would also be interesting to evaluate the 
behavior of the routing algorithm to minimize the cycle time, the 
work-in-process (WIP) or other performance measures. Topics for 
future research on the area include the routing in queueing 
networks with cycles, e.g., to model many important industrial 
systems that have reverse streams of products due to re-work, or 
even the extension to GI/G/c/K queueing networks. 
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